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Current methods of evolving a spacetime containing one or more black holes are plagued by in- 
stabilities that prohibit long-term evolution. Some of these instabilities may be due to the numerical 
method used, traditionally finite differencing. In this paper, we explore the use of a pseudospectral 
collocation (PSC) method for the evolution of a spherically symmetric black hole spacetime in one 
dimension using a hyperbolic formulation of Einstein's equations. We demonstrate that our PSC 
method is able to evolve a spherically symmetric black hole spacetime forever without enforcing 
constraints, even if we add dynamics via a Klein-Gordon scalar field. We find that, in contrast to 
finite-difi'erencing methods, black hole excision is a trivial operation using PSC applied to a hyper- 
bolic formulation of Einstein's equations. We discuss the extension of this method to three spatial 
dimensions. 

04.25. Dm, 02.70.Hm 



in 



> ■ 

vn 
o 
in 

o : 
o . 
o . 

"o '. 

I 



X 



I. INTRODUCTION AND SUMMARY 

A major thrust of research in classical general relativ- 
ity in the past decade has been to devise algorithms to 
solve Einstein's equations numerically. Despite advances 
in our analytic understanding of general relativity, we 
still do not know what all the features of the theory re- 
ally are. Numerical solutions will continue to provide 
fresh insights into the theory as they have in the past, 
for example, critical behavior in black hole formation |^] 
and the formation of toroidal black holes ||] . 

New urgency has been injected into numerical rela- 
tivity by the imminent deployment of LIGO. The prime 
target for LIGO is coalescence of binary neutron star and 
black hole systems. The waveform is reasonably well pre- 
dicted by the post-Newtonian approximation when the 
binary components are at large separation. However, ex- 
tracting the most important physics requires us to be able 
to deal with fully non-linear general relativity as the sys- 
tem spirals together and coalesces. Moreover, a number 
of people believe that there is a significant event rate for 
the coalescence of massive black hole systems (~ 20Mq) 
In this case, LIGO is most sensitive to waves emitted 
from the strong field regime. Indeed, without some the- 
oretical guidance as to what to expect from this regime, 
it is possible we may miss these events entirely Q . 

However, the goal of developing a general algorithm 
that can solve Einstein's equations for two black holes 
has remained elusive. All attempts to date have been 
plagued by instabilities. These instabilities are caused by 
an interplay of three factors: (1) Einstein's equations are 
an overdetermined system, with the evolution equations 
subject to constraints. So if, for example, you choose to 
solve only the evolution equations, then there can be un- 
stable solutions that are in fact solutions of the evolution 
equations, but do not satisfy the constraints. Small nu- 



merical errors may cause these solutions to appear and 
swamp the true solution ( "constraint- violating modes"). 
(2) The coordinate freedom inherent in the theory means 
that it is very easy to impose coordinate conditions that 
lead to numerical instabilities ("gauge modes"). (3) Ex- 
perience has shown that the kind of boundary conditions 
we choose and how we implement them can affect the 
stability of an algorithm enormously. 

Similar instabilities have hampered efforts to solve the 
related problem of binary neutron stars; only very re- 
cently has there been some success in finding stable 
algorithms. However, black hole evolutions face an addi- 
tional obstacle that is absent in the case of neutron stars: 
for neutron stars the gravitational field is everywhere reg- 
ular, but for black holes one must somehow deal with the 
physical singularity that lurks inside each hole. 

There are two main approaches for handling these sin- 
gularities. The first is to use gauge conditions (e.g., max- 
imal slicing) that avoid the singularities altogether. Such 
conditions, however, lead to large gradients in the grav- 
itational field variables near the horizon. These grow 
exponentially in time and ultimately cannot be resolved 
by the numerical evolution, causing the code to crash. 
The alternative approach is to excise the region contain- 
ing the singularity from the computational domain and 
evolve only the exterior region. If the excision boundary 
is placed inside the horizon of the black hole, causal- 
ity assures us that we do not need to impose a physical 
boundary condition there. 

However, black hole excision is only known to be math- 
ematically well-posed if the evolution equations are hy- 
perbolic with characteristic speeds less than or equal to 
c. In this case, the structure of the equations guaran- 
tees that even unphysical modes present in the solu- 
tion (gauge modes, constraint-violating modes) behave 
causally and cannot propagate out of the horizon. For 
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many representations of general relativity such as the 
usual ADM formulation, the evolution equations are 
of no mathematical type for which well-posedness has 
been proven, so the suitability of these formulations for 
black hole excision must be determined empirically on 
a case-by-case basis. It is in part for this reason that 
much attention has been recently focused on hyperbolic 
representations of Einstein's equations ]^-[l9|] . 

It is still unclear whether hyperbolic formulations are 
computationally advantageous. However, it has been 
shown that the formulation of the evolution equations 
can affect stability. For example, some instabilities can 
be eliminated by changing from one formulation of Ein- 
stein's equations to another ||2^] or by modifying the evo- 
lution equations to change the spectrum of unphysical 
modes |2l). 

Nevertheless, it is likely that many instabilities encoun- 
tered in practice are due to the numerical implementation 
of the evolution equations and of the boundary condi- 
tions. Even for well-understood systems, it is far from 
guaranteed that any given numerical approximation to 
the continuum equations will be stable. The well-known 
Courant instability is a trivial example. Hence it is pru- 
dent to explore alternative numerical methods that may 
offer a shorter path to the goal of long-term stability. 

Traditionally, black hole spacetimes have been evolved 
using finite-difference (FD) methods. Current FD codes 
for evolving black hole spacetimes with excised horizons 
are mostly based on a numerical technique known as 
causal differencing |p^-p7|, which allows one to update 
the fundamental variables in time while avoiding numer- 
ical problems associated with superluminal grid speeds. 
This technique has been used successfully to propagate 
an excised hole across a grid, even when grid points fall 
into or emerge from the horizon p5| . 

However, causal differencing is complicated because it 
requires interpolation, and it has to deal with points 
"missing" in some irregular fashion near the excision 
boundary. The FD operator that performs the interpo- 
lation depends on the shape of the excision boundary. 
Furthermore, even for a given excision boundary and a 
given target point, the operator is not unique. Hence 
one must construct a large number of interpolation op- 
erators, each of which involves some arbitrary choice, in 
order to perform interpolations on the entire grid. It is 
only by trial and error that one finds operators that result 
in a stable evolution scheme. For simulations of a sin- 
gle spherical black hole on a Cartesian three-dimensional 
grid, we have found a case in which changing a single in- 
terpolation operator that is used only for a single target 
point on the grid makes the difference between a stable 
and an unstable code. 

Another limitation of FD methods is the difficulty of 
imposing boundary conditions at the outer boundary of 
the calculation. There are two aspects to this problem: 
First, one must formulate a procedure for handling the 
boundary, such as imposing an analytic condition (e.g., 
a Sommerfeld condition) on the fundamental variables. 



matching to a wave perturbation described by the Zer- 
illi equation |^^, or matching to a characteristic evolu- 
tion code that propagates the solution out to null infinity 
p9|-^. Second, one must construct a FD approximation 
of cither the analytic boundary condition or the matching 
condition. It can be difhcult to find such an approxima- 
tion that yields a stable evolution. 

In this paper we explore an alternative computational 
strategy: a pseudospectral collocation (PSC) scheme. In 
our PSC evolution scheme, the solutions to a set of hy- 
perbolic differential equations are approximated as series 
expansions in a set of orthogonal basis functions (e.g., 
Chebyshev polynomials) in space, and these coefficients 
are integrated forward in time using the method of lines. 
PSC has three important advantages over FD: (1) No 
ad hoc interpolation operators are required to determine 
field values at an arbitrary point because the solution 
provided by PSC is an analytic function given every- 
where on the computational domain. (2) Boundary con- 
ditions are imposed directly on the basis functions, with 
no approximations, in a straightforward manner. (3) For 
smooth solutions, PSC will converge to the actual so- 
lution exponentially as the number of basis functions is 
increased. A FD solution, on the other hand, never con- 
verges faster than algebraically with the number of grid 
points. Thus, for a given accuracy, PSC requires far less 
CPU time and memory than FD methods. 

PSC has been applied successfully to solve problems in 
many fields, including fluid dynamics, meteor olog y, seis- 
mology, and relativistic astrophysics (cf. p^-pSl). For 
example, PSC has been applied successfully to model 
stellar core collapse ||3^ ] and construct equilibrium se- 
quences of irrotational binary neutron stars [ p7| . For 
black hole spacetimes PSC has been applied successfully 
to solve initial data for the standard field equations |3^] 
and the conformal field equations |^^, to find apparent 
horizons [ [40[, to solve the shift vector equation for a Kerr 
black hole ||4l| , and to evolve Einstein's equations in null 
quasi-spherical coordinates ||4^ . 

Here we evolve a spherically symmetric black hole 
spacetime in one spatial dimension by applying PSC 
methods to a hyperbolic formulation of Einstein's equa- 
tions, the "Einstein-Christoffel" (EC) system A hy- 
perbolic formulation provides a well-defined prescription 
for imposing boundary conditions. At a boundary, the 
fundamental fields can be decomposed into characteris- 
tic fields, which in the case of the EC system propagate 
either along the light cone or normal to the spatial folia- 
tion. Boundary conditions are imposed on the incoming 
(with respect to the computational domain) characteris- 
tic fields, but not on the outgoing fields. After imposing 
the boundary condition, the fundamental fields are re- 
constituted from the characteristic decomposition. If the 
excision boundary is placed inside the event horizon of 
a black hole, then (for appropriate coordinate systems) 
all characteristic fields are outgoing at this boundary, so 
no boundary conditions are needed there. Thus, black 
hole excision is a trivial operation. At the outer bound- 
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ary of the domain, boundary conditions corresponding to 
no incoming radiation can be imposed on the incoming 
characteristic fields. 

We find that our PSC method is able to evolve a spher- 
ically symmetric black hole spacetime forever. Further- 
more, we find that the solution converges exponentially to 
the exact solution as the number of basis functions is in- 
creased. We discuss the time-stepping algorithms, outer 
boundary conditions, and gauge conditions required for 
stable evolution and how these depend on the particular 
slicing of the Schwarzschild geometry we wish to repro- 
duce. We also show that our PSC method can handle 
dynamics by evolving a black hole spacetime containing 
a scalar field. Finally, we outline a strategy for applying 
this method to two black holes in three spatial dimen- 
sions using multiple domains, and we present tests of 
this domain decomposition idea in spherical symmetry. 

In Sec. |l| we list the basic equations we use to evolve 
a spherically symmetric spacetime, boundary conditions, 
gauge conditions, initial data, and diagnostics. In Sec. 
Ill we introduce PSC and describe our numerical meth- 
ods. In Sec. IV we present numerical evolutions of single 
black hole spacetimes in spherical symmetry. Finally, in 
Sec. we discuss our results and the generalization of 
our methods to multiple dimensions and to binary black 
hole spacetimes. 



II. BASIC EQUATIONS 

A. Einstein-Christoffel system 

We adopt the "Einstein-Christoffel" (EC) hyperbolic 
representation of Einstein's equations M] . In this formu- 
lation the evolution equations are written in first-order 



symmetric hyperbolic form, and all characteristic curves 
are directed either along the light cone or normal to the 
spatial foliation. The EC system takes as fundamental 
quantities the familiar three-metric and extrinsic curva- 
ture plus only 18 additional "connection" variables (in 
three spatial dimensions). Unlike some hyperbolic repre- 
sentations of general relativity, the EC system requires 
no derivatives of the stress-energy tensor. 
We write the metric in the usual 3+1 form 



-N^dt^ + gij{dx' + I3'dt){dx^ + (i^dt), (2.1) 



where gij is the three-metric, /3* is the shift vector, and 
N is the lapse function. In the EC formulation it is not 
the lapse function N that is freely specifiable; instead, 
one arbitrarily prescribes the densitized lapse function 
a, defined by 



N 



(2.2) 



where g is the determinant of the three-metric. The use 
of a densitized lapse does not fix the temporal gauge free- 
dom in any way: one can in principle obtain any lapse 
function N by an appropriate choice of a. 

To write down the EC evolution equations, first define 
the new variables 



(2.3) 



where r\j is the affine connection associated with gij, 
and parentheses and brackets denote symmetrization and 
antisymmetrization, respectively. The quantities fkij will 
be taken as fundamental variables along with and the 
extrinsic curvature Kij. 

The EC evolution equations can be written in the form 



dog.j - -2NK 



doL 



kij 



Nd.K, 



N{g*'\KuKij - 2KkiKij) + 5*"''g""(4/fo„i/[;„]j + ■ifkm[nfl\ij - fikrnfjln 
+ ^f{ij)kf[ln\m + '^fkm{ifj)ln — ^fklifmnj + 20fkl(ifj)mn ~ 13/ifcz/jmn) 

- d^dj In a - {di In a) {dj In a) + 25^ g'^'g™" {fkmndi In a - fkmidn In a) 
+ 9''''[{'^f{ij)k ~ fktj)di Ina + AfkK^^dj) \na~ 3{f,kidj lna + f^kid. In a)] 

- SirSij + AngijT}, 

N {g"^'^[AKk(ifj)mn - '^fmn(iKj)k + Kij{2fmnk " 3/fem„)] 

-f 2g""gP9[i^™p(5fc(J,),„-2/, 

n(i9j)k) ~^ 9k{i-^j)rn(.^fnpq ^fpqn) 
+ Kmn{'^fpq(i9j)k " ^9k(if])pq)] ~ K^jdk In Q; 

+ 2g'^'^{K^(^^gj)kdn Ina - Kmngk(idj) In a) + lQT:gk(iJj)}- 
Here the symbol 9o is the time derivative operator normal to the spatial foliation, defined by 

do = dt- £i3, 

where £ denotes a Lie derivative. The matter terms are 



(2.4a) 



(2.4b) 



(2.4c) 



(2.5) 
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T = 




(2.6b) 


J^ = 


ll ^fil^J 


(2.6c) 


Sij = 




(2.6d) 



where T^^ is the stress-energy tensor, ^'^^ g^u is the four- 
metric, 71^ is the unit normal to the spatial foliation, and 
7j'' is the spatial projection operator n'^n^ -I- '^^''g/ ■ 

Not e th at for each pair, the evolution equa- 

tions (2.4) can be written in the form (dropping the i 
and j indices on gij , Kij , and fkij ) 



dag ■ 



-2NK, 



(2.7a) 



dofk + NdkK = S, 



(2.7b) 
(2.7c) 



where R and 5 are nonlinear terms that contain no 
derivatives of the fundamental variables (but may con- 
tain spatial derivatives of the arbitrary gauge function 
a). Except for the right-hand sides, equations (2.7a) are 
just Og = written in first order form. Thus one can 



think of the EC system (2.4) as a set of six (one for 



each {i,j} pair) coupled quasilinear scalar wave equa- 
tions with nonlinear source terms. 



A solution of the evolution equations ( 2.4) is not a solu- 
tion to Einstein's equations unless twenty-two constraints 
are also satisfied. These are the Hamiltonian constraint 



+ /^jfe(8/™/„ - 20/,„„)]} + IQiip = 0, (2.8a) 
the three momentum constraints 

C, = 5'''{ff""[-^.:fc(3/;mn - 2/™„/) - Kk^hn] + d,Ku - duKa} + Sir J, = 0, (2.8b) 
and the eighteen constraints 

Ckij = dkgij - 2fkij + ^g^"^{fim{ig])k - gk(ifj)im) = (2.8c) 



that relate fkij to spatial derivatives of the three-metric. 
If the constraints are satisfied for the initial data, they are 
preserved by the evolution equations for all time. As for 
any formulation of Einstein's equations, however, numer- 
ical approximations may spoil this constraint-preserving 
property. 



B. Reduction to spherical symmetry 

The most general spherically symmetric metric can be 
written in the form 



ds^ = -N^dt^ + grr{dr + (3'' dtf 
+gTr^ {de^ +sin^ 0d(j>'^), 



(2.9) 



where the transverse metric component is defined by 

_ goe _ g^4> 



sm 6 



(2.10) 



The two nonvanishing independent components of the 
extrinsic curvature are K^r and the transverse extrinsic 
curvature 



Kt 



Kee Kit, 



sin^ I 



(2.11) 



The nonvanishing components of fkij are frrr, the trans- 
verse component 



frl 



free feer fr4 



J.2 r2 sin^ 6 2r'^ sin^ I 

and the additional components 

frr9 — grr COt 6, 

feee = 2r^5T cot 9, 

- r^gT sin ^ cos 0, 
: r^gr sin 6* cos 6*. 



(2.12) 



(2.13a) 
(2.13b) 
(2.13c) 
(2.13d) 



The e voluti on equations for these additional compo- 



nents (2.13) are automatically obeyed if the evolution 
equations for the metric are satisfied. We therefore do 
not treat these quantities as independent variables, and 
wherever they appear in the equations we replace them 
with t he a ppropriate metric components using equa- 
tions ( 2.13 ). Of course, all angular dependence due to 
these terms drops out. 

We therefore take as fundamental variables the six 
quantities grr, gr, Krr, Kt, frrr, and frT- Using (pl|), 
we obtain the following evolution equations for these vari- 
ables: 
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2NKrr + 2grrdrP'', 



dt9T - P^'drgT = ~2NKt + 2'—gT 

r 



N 

dtKrr — P^drKrr H drfr 

grr 



N 



N 



rr 



1 4/, 



gr 



6 



fr 



{dr Ina)^ + ( f^^^ ] drill a 



dtKr - (3'drKT + —drfrT = N ( KtK\ + 
grr 



1 frT 



gr J \gT 

+ 2Krrdrl3'' + inNiTgrr - 2Srr) 

2(3' 



grrgr grr 

r TCr^ f 

dtfrrr — P^drfrrr + NdrKrr = N 

+ ifrrrdrP'' + ffr-r^^/?'' + 167riV J, 

frT 



dr In (5 



(2.14a) 
(2.14b) 

9^ In 5 

(2.14c) 
(2.14d) 



4g„— f 3^ - + l-dr\na\- Krr f 10^ + frr ^ - + 9. Ill 5 



gr \ gr r 
^2 ar 



dtfrT - P'' drfrT + N BrKr = iV 



Kt 2 



gT 



f^rr ~ 1^ ^ 



gT 



drP^ + — ) /.T. 

r 



Here 



N 



a = ar sm ( 



gTJgrr 



and we have explicitly included the terms involving the Lie derivative of the shift vector. 
The six fundamental variables obey four constraints that can be obtained from (p.8[): 



C 

c 

CrT 



_ drfrl 



1 



frT / 2 7f, 



rT 



grrgT '2r^gT grrgT \r 2gT 
drKr , 2Kt frT ( ,^r , 

\ A ^ H 

gT rgx gT V 9T 



-f\ 
AnJr = 



gr V 25t/ 



drgr 



grr frT 

gT 



^frrr — 0; 



drgT + ^ - 2frT = 0. 

r 



(2.14e) 
(2.14f) 

(2.15) 



(2.16a) 
(2.16b) 
(2.16c) 
(2.16d) 



We do not explicitly solve the constraints during our evo- 
lution, but instead we use them as error estimators. 



C. Boundary conditions 

Boundary conditions are imposed on the above evo- 
lution equations (2.14) via characteristic decomposition. 
Consider a first-order symmetrizable hyperbolic system 



dtU + A'd,U = R, 



(2.17) 



where U is the vector of variables, i? is a vector, and the 
three A* are matrices. Then for a particular unit vector 
^i, the solutions Uc to the eigenvalue problem 



(2.18) 



define the characteristic fields normal to the direction , 
and the eigenvalues Vc define the characteristic speeds of 
these fields. Each of the characteristic fields Ur can be 



thought of as a plane wave solution moving in the direc- 
tion with speed Vc- One is allowed to impose boundary 
conditions only on characteristic fields that propagate 
into the computational domain, but not on fields that 
propagate out of the domain. 



For the evolution equations (2.14), the characteristic 
fields in the radial direction (^,. = ^/g^) are 



(2.19a) 
(2.19b) 





— grr 




-n, 




= gT 




-n, 




frrr 


(Wc = 


-f3-±agT), 


= Kt± 


frT 
■\J Qrr 


(Vc = 


-p'^±agT). 



(2.19c) 



The characteristic speeds of the metric variables corre- 
spond to propagation along the timelike normal to the 
foliation and the characteristic speeds of the other quan- 
tities correspond to propagation along the light cone. 
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Thus, if the inner boundary of our domain moves along a 
spacehke trajectory, all characteristic speeds are negative 
(with respect to r) there, so no boundary conditions need 
to be imposed. At the outer boundary, boundary condi- 
tions are imposed only on those quantities with negative 
characteristic speeds (usually C/°, [7°, , and U^). 

We have experimented with three types of outer 
boundary conditions. The first, which we call the freezing 
boundary condition, is 



dtUc = 



(2.20) 



applied to all incoming characteristic fields Uc- This cor- 
responds to no incoming radiation at the boundary; how- 
ever, for nonlinear or inhomogeneous problems this is not 
strictly correct unless the boundary is at infinity. 
The second is the Robin condition 



drK{U-Uoo)]=0, 

which assumes that U behaves like 

const 



(2.21) 



(2.22) 



at large r. For a given incoming variable Uc, appropriate 
values of the parameters Uoc and n can be found from the 
analytic representations of the Schwarzschild geometry in 
section II D below. 

Finally, the constraints ( 2.16| ) can be used to derive 



mixed Neumann-Dirichlet boundary conditions f or fou r 
of the characteristic fields: Equations ( 2.16c ) and ( 2.16d ) 
can be used directl y as bo undar y con ditions on [/^ and 
U^, and Equations ( ^.16^ ) and ( ^.16bD can be combined 
to yield boundary conditions on and U.^ . We use only 
three of these boundary conditions — the ones for U^, [7°, 
and — because is outgoing at the outer bound- 
ary and therefore needs no boundary condition there. 
Similarly, in three spatial dimensions the constraints can 
be used to derive twenty-two relations among the thirty 
characteristic fields, some of which can be used as bound- 
ary conditions on the incoming fields. 

We describe our numerical implementation of bound- 
ary conditions in section III D , 



D. Coordinate systems 



It is convenient to choose a coordinate system in which 
the Schwarzschild geometry is time-independent. Fur- 
thermore, since we wish to include the apparent horizon 
in our computational domain, we must choose coordi- 
nates such that the spacelike slices labeled by constant 
values of coordinate t penetrate the horizon and are non- 
singular there. Here we list several coordinate systems 
that satisfy these properties. 



1. Kerr-Schild coordinates 



In this coordinate system, also referred to as ingoing 
Eddington-Finkelstein coordinates [^3|| , ingoing null rays 
have unit coordinate speed. In addition, the radial coor- 
dinate r is chosen such that 47rr^ is the surface area of 
a sphere at that radius. In this coordinate system the 
Schwarzschild solution takes the form 



a = 



9rr 1 

9T = 1, 

H 

Kt 

frrr 
frT 



2M 



2M 

r 

1 + 



2M 
r 



2M 

2~ 



2M 

i|4. 

r 

1 

7 

r 



i + H 

r 

2M 
r 

7M 
r 



1 - 

1/2 



-1/2 



(2.23a) 
(2.23b) 

(2.23c) 
(2.23d) 
(2.23e) 
(2.23f) 
(2.23g) 
(2.23h) 



where M is the mass of the hole. The event horizon is 
coincident with the apparent horizon and is located at 
r = 2M. 



2. Painleve-Gullstrand coordinates 



In this coordinate system [E6 44-4q| the spatial three- 



metric is flat and the Schwarzschild solution is particu- 
larly simple: 



9rr 

9t 
a 

13'' 
Kt 

frrr 
frT 




2M 



(2.24a) 
(2.24b) 
(2.24c) 

(2.24d) 

(2.24e) 

(2.24f) 
(2.24g) 
(2.24h) 



The horizon is again located at r = 2M. 
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3. Harmonic time slicing, areal radial coordinates 

If one requires the time coordinate to satisfy Ut = 0, 
the radial coordinate r to correspond to the areal radius, 
and the coordinate system to be regular at the horizon, 
then the Schwarzschild solution takes the form 

4Af2 



9rr 
9T 



= 1 



1, 



5=1 



2M 
r 



4M 



-f2 



Kt 

frrr 
frT 



4Af2 r 



5 2 + 



3M 4M2 4M3\ 



4 7Af 12A-/2 
1 



20Af^ 



(2.25a) 
(2.25b) 

(2.25c) 
(2.25d) 
(2.25e) 

(2.25f) 
(2.25g) 
(2.25h) 



The horizon is at r = 2Af. 



4- Fully harmonic coordinates 

The Schwarzschild solution can also be written in a 
coordinate system where all coordinates satisfy Dx^ = 
and are regular at the event horizon M,E8|: 



grr — \ + e + + e", 



9t 



M 
r 

M 



M 



1 + ^ (1 + ^^)1 + ^ 



3M 
r 

3M 

r 



(2.26a) 
(2.26b) 

(2.26c) 

(2.26d) 



Krr = 

Kt = 

frrr ^ 
frT = 

where 



9T V 2 

4A//2 ^ 



,,3 



1 



1 M 
- + 72 



1(1 + a 



-^(ll-2e + 962) 



2A/ / M 
1 + — 

r \ r 



Here the horizon is located dX r — M . 



(2.26e) 

(2.26f) 
(2.26g) 
(2.26h) 

(2.26i) 



E. Einstein-Klein-Gordon system 

To add dynamics to the spherically symmetric prob- 
lem, we introduce a Klein-Gordon scalar field (j) with 
stress-energy 



Defining the quantities 



n = -N-^do^, 



the matter terms (|2.6| ) are given by 



(2.27) 



(2.28a) 
(2.28b) 



(2.29a) 

(2.29b) 
(2.29c) 

(2.29d) 



The scalar field obeys □</> — 0, which in spherical sym- 
metry takes the form 



n 



\ V dr\na 

9rr gr J grr V 5T T 



dt^r - P'^dr^r + Ndrli = - NI{ I — - - - -f 9, In 5 ) + ^rdrP"" 

grr gr r 



(2.30a) 
(2.30b) 



Only the derivatives of the scalar field <p appear in equations (2.27) and (2.3C), so ip itself need not be evolved. 
For the evolution equations (2.30), the characteristic fields in the radial direction (^^ = \fgrr) are 



f/^ = n± 



(2.31) 



with characteristic speeds — /?'" ± figT- 
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F. Apparent horizon 



1. Algebraic conditions 



A marginal outer trapped 2-surface is defined by the 
equation 



9- 



0, 



(2.32) 



where Di is the covariant derivative compatible with the 
three-metric, and s* is the outward-pointing spatial unit 
normal of the surface. In spherical symmetry, equa- 



tion (2.32) reduces to 

frT 



/gr 



0. 



(2.33) 



The apparent horizon is the outermost surface at which 
( 2.33 ) is satisfied. On each time slice, the co ordin ate ra- 
dius of the horizon rah is located by solving (2.33) using 
a standard root-finding algorithm. 

If the horizon is to remain at a fixed coordinate radius 
as the spacetime evolves, the following relation must be 
obeyed at the horizon: 



JJ^ _ 1 + Sirr'^gTip - Jr/y/gT^) 
agr 1 - Snr^grip - Jr/ ' 



(2.34) 



Equati on ( 2.34) can be derived by setting the time deriva- 
tive of ( 2.33| ) equal to zero and substituting the evolution 
and constraint equations to eliminate time and spatial 
derivatives. 



G. Gauge conditions 

The question of which gauge conditions one should im- 
pose on a numerically generated spacetime is one of the 
key unsolved problems in numerical relativity. In princi- 
ple, the coordinate invariance of general relativity allows 
one to make this choice arbitrarily. However, a poor 
choice may not only obscure the physics one is searching 
for in the simulation, but may also allow rapidly-growing 
gauge modes that halt the code altogether. 



_ Krr9T 
v/ grr 



2/, 



rT 



9T 



dra + ay/g^gr 



drfrrr 
grr 

drKr 



The simplest gauge choices we consider here are alge- 
braic ones: we set a and /3* equal to their analytic values 
for some parameterization of the Schwarzschild solution 
in section II D| that we wish to reproduce numerically. 
While these gauge conditions are obviously applicable 
only to test problems, they provide a simplified setting 
in which to study the properties of our evolution scheme. 
Such conditions have been used extensively in 3D test 
problems [p5| . 

Next we consider other algebraic gauge conditions that 
are independent of a particular analytic solution, but are 
still of limited generality. For instance, one might re- 
quire that the radial coordinate remains areal, or in other 
words, that gx is time- independent. Using our variables, 
this condition can be written 



13'' frT - agry/g^KT = 0. 



(2.35) 



One might also require that the ingoing coordinate speed 
of light takes on a prescribed value c_ : 



agT + /3'' + c_ = 0. 



(2.36) 



These conditions are not generalizable to two black holes 
because the first relies on the notion of an areal radial co- 
ordinate and the second assumes that a unique "ingoing" 
direction exists at every point in spacetime. Nevertheless, 
these and similar gauge conditions have proven useful for 
studies of single-black-hole spacetimes [i3p9t|. Further- 
more, if imposed only at one point, they can be used 
as boundary conditions on more general elliptic gauge 
choices, described below. 



2. Elliptic conditions 



We also explore gauge conditions that should be ap- 
plicable to general spacetimes. For the shift vector, we 
consider two elliptic equations due to feOj : minimal strain 



idrfrT _ 2 



2/rT / 5/r 



9r 



gT V 9rr J 9T \ 9rr 

'^KxfrT ^ Krr f frrr ^ 2 ^ &frT 
{9tY 9rr V 9rr T gT 



frT_ 

9t 

= 0, 



(2.37) 



which minimizes changes in the three-metric in a global sense, and minimal distortion. 



ff^pr^ j frrr 



2/, 



rT 



9t 



dj-frrr ^^r fi 



rT 



9T 



+ 



Kt 
9t 



Kr 
9ri 



Qrr 



drKr 
grr 



frrr \ 


^ frT 




5frT 10 


grr / 


9T 


V 9rr 


9T r 




( frrr 




Krr f frrr ^ 


9t 


\ grr 




grr \ grr 



9t 



0, (2.38) 



which similarly minimizes changes of the conformal three-metric. 
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For the densitized lapse function, we consider the stationary mean curvature condition dtK = 0, which in terms of 
our variables can be written 



frrr 
Qrr 

9t 



2fr 



9r 



4 

9T r 

Krr 



dra 4 

9T 



rgr 



9T 



frrr 
9rr 



frT 

9t 



frrr 
9rr 



5/rT 
9T 



6 



drKr 



2drKT ^ AKt (I frT 



9t 



9t 



9t 



AngrrT 

^ 2K„ (4. frT 
9rr V 9T 



frn^ 



= 0. 



(2.39) 



This condition was also discussed by [|0[, and is best 
known in the special case = 0, when it reduces to the 
familiar maximal slicing condition. Use of the stationary 
mean curvature lapse combined with either the minimal 
strain or minimal distortion shift vectors has been re- 
cently encouraged by ||5l|,^. 

Each of the elliptic equations ( |2.37D , ( |2.38| ), and ( |2.39| ) 



requires two boundary conditions. At the horizon, we 
either set a and /?'' t o prescribed values, or we im- 
pose ( 2.34 ) and ( 2.36 ). At the outer boundary, we 



again can set a and (3^ to prescribed values, we can im- 
pose ( 2.35 ) a nd (2.36 ), or we can impose Robin conditions 
of the form ( 2.21 ) on a and /?''. 



H. Mass 

It is useful for diagnostic purposes to compute the mass 
of the spacetime. In spherical symmetry, the total mass 
inside an invariant spherical surface labeled by coordinate 
r is well-defined and given by the Misner- Sharp formula 
[B3| , which for our variables reads 



MMs(r) 



r\/9T 



1 



9T 



rT 



(2.40) 



III. PSEUDOSPECTRAL COLLOCATION 
METHODS 

A. Introduction 

Consider a system of L evolution equations of the form 

(3.1) 

for 1 < £ < L, where {f^^\x,t)} is the solution, and 
j:"!^) [{ are (possibly nonlinear) functions of {/'■^-'} 
and their spatial derivatives. Approximate each func- 
tion /(^^ of the solution as a finite sum of basis functions 



f^;^{x,t) = Y^fi\t)4>f{x). 



(3.2) 



fc=0 



For smooth functions as ^ cx) the approximation is 
exact. Corresponding to the approximate solution {Z^"*} 
is a residual 



(^) 

N 



IN 



(3.3) 



for each evolution equation. 

In PSC the spectral coefficients fl^\t) are determined 
by demanding that the residuals vanish at a fixed 
set of N collocation points Xn. In other words, it is de- 
manded that the system of differential equations (3.1) 
be satisfied exactly at the collocation points {af„}. The 
choice of the collocation points is intimately related to 
the choice of basis functions used in the approximate so- 
lution. In the following subsection, we discuss how they 
are chosen. 



B. Expansion basis and collocation points 

For the remainder of this section, we will restrict our- 
selves to problems with one spatial dimension (ID). The 
choice of an expansion basis depends upon the particu- 
lar problem being solved. For example, the natural ex- 
pansion basis for a ID problem with periodic boundary 
conditions is a Fourier series. For more general bound- 
ary conditions, such as the ones we will impose in our 
black hole evolutions, Chebyshev polynomials are a ro- 
bust choice for the basis functions. Chebyshev polyno- 
mials are defined on the interval 



by 



I =[-1,1] 



Tfc(a;) = cos(fccos ^ x). 



(3.4) 



(3.5) 
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A function / on I is approximated as |^ 

N 

/jv(a.,i) = ^A(i)Tfe(a;). (3.6) 

fc=0 

Note that in order to use this expansion, we must specify 
a mapping from our physical domain [rminj^'max] to I. 
The simplest choice is a linear mapping, but other choices 
may work better. 

For a Chebyshev expansion, a convenient choice of the 
collocation points is 



Xr, — COS ■ 



N 



(3.7) 



At these collocation points, the Chebyshev polynomials 
satisfy the discrete orthogonality relation 



2 ^ 1 

5jk = -ITj— —Tj{Xn)TkiXn), 



where 



Ck 



2, k = Oor N 
1, 1 < fc < iV- 1. 



(3.8) 



(3.9) 



Using the orthogonality relation, the spectral coefficients 
are given by 



Nck 



N 

E, 

n=0 



1 



-fN{Xn)Tk{Xn)- 



Since 



Tk{Xn) 



COS ■ 



irkn 



(3.10) 



(3.11) 



fast cosine transforms can be used to comp ute ( p^ ) at 
the collocation points and to evaluate ( |3.1C| ). 

In PSC, the focus is not on the set of spectral coeffi- 
cients {/fc(i)}, but on the equivalent set {/(x„,f)}, the 
approximate solution evaluated at the collocation points. 
In particular, the approximate solution to (3.1) would be 
given by evolving 



dtf^j^\xn,t) ^ T'^'\xn,t), 



(3.12) 



for 1 < £ < L and < n < A^. Given initial con- 
ditions f^ ^\x^ ), and appropriate boundary conditions. 
Equation (3.12) can be evolved forward in tim e using the 
method of lines, described in section HID. Since the 



focus is on grid-point values, and not the spectral co- 
efficients, it is possible to reuse large amounts of code 
developed for FD methods. 



'^For Chebyshev bases the conventional notation is that k 
runs from to N , not — 1; thus, there are A'^ -I- 1 coefficients 
and collocation points. 



C. Computation of derivatives 

T he m ain differences between PSC and FD in evolv- 
ing (3.12) are the choice of collocation (grid) points Xm 



how spatial derivatives are computed, and how boundary 
conditions are imposed. In PSC, spatial derivatives are 
computed analytically from the series expansion 



dfN{x,t) 
dx 



k=0 



dx 



(3.13) 



This derivative can be written as another sum over 
Chebyshev polynomials 



dfN{x,t) 
dx 



N 



by using the simple recursion relation 

Cfc/fc(i)-/fe+2W + 2(fc + l)/fe+i(t), 

where 



Ck = 



fc = 
fc > 1. 



(3.14) 



(3.15) 



(3.16) 



Evaluating a derivative requires two fast transforms; 
the first to comput e the spectral coefficients neede d in t he 
recursion relation (3.15), the second to evaluate (3.14). 



D. Time evolution and application of boundary 
conditions 



We evolve our hyperbolic system using the method 
of lin es. I n this method, we cast our system into the 
form ( ^.12| ) and use a standard ODE solver to integrate 
the equation in time. For most of the results presented in 
this paper, we have used a fourth-order explicit Runge- 
Kutta method. One of the drawbacks of using PSC is 
that the Chebyshev collocation points (^^) are clustered 
near the domain boundaries. This places a more severe 
Courant stability limit At ^ OiN-"^) on a wave equa- 
tion than for FD, where At ~ Ax ~ 0{N~^). Because 
of the superior spatial convergence of PSC, however, this 
restriction is not as severe as it may seem at first glance. 
In fact, to retain the accuracy gained by the spatial res- 
olution, it may be necessary to use a time step smaller 
than that demanded by stability. Moreover, if the sta- 
bility restriction on the time limit becomes too severe 
for practical evolutions, one can implement implicit or 
semi-implicit time-stepping schemes. 

One of the advantages of PSC over FD is in how 
the boundary conditions are applied. In FD, derivatives 
are approximated by differences of field variables at grid 
points. The pattern of grid points used must typically be 
modified at the boundaries of the numerical grid. Con- 
sequently, boundary conditions can be difficult to formu- 
late in FD. In PSC, on the other hand, the approximate 
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solution is given over the entire domain. As seen in the 
previous section, derivatives are computed analytically; 
therefore nothing special needs to be done to compute 
the derivative at a boundary. Furthermore, since there 
are collocation points on the boundary of the domain, the 
application of boundary conditions is straightforward in 
PSC. One simply demands that the approximate solution 
satisfy the exact boundary condition at the boundary col- 
location point. 

The boundary conditions are applied during the time 
step by modifying .F*^^-* [{/*^^^}] (cf. equation 3.1) at the 
boundary points so that the boundary conditions are sat- 
isfied. In this paper we are interested in applying bound- 
ary conditions on a hyperbolic system of evolution equa- 
tions. As described in section [IC, the solution to a 



hyperbolic system can be written in terms of characteris- 
tic fields that propagate with corresponding characteris- 
tic speeds. Physically we know that boundary conditions 
need only be applied to the incoming characteristic fields. 

Therefore, to impose a boundary condition at a domain 
boundary x = Xb, we first compute the time derivatives 
of the characteristic fields Uc{xb) at the boundary. We 
then apply a boundary condition to those fields that are 
propagating into the domain; the remaining fields are un- 
touched. Finally, we reconstruct the time derivatives of 
the fundamental variables at Xb and use these values in 
the time update. Failure to impose a boundary condi- 
tion on an incoming field or imposition of a boundary 
condition on an outgoing field almost always leads to an 
unstable evolution. 

A Dirichlet boundary condition Uc{xb,t) — g{t) is ap- 
plied by ensuring that the time derivative of the incoming 
characteristic field at the boundary collocation point is 
set to dg/dt. A boundary condition such as Neumann or 
Robin that involves the spatial derivative of the charac- 
teristic field is enforced by replacing the spatial derivative 
at the boundary with the appropriate value in order to 
satisfy the boundary condition. 



p4| to handle communication between multiple domains 
and for parallelization of our code. 

The extension of our method from one domain to mul- 
tiple domains is straightforward. We evolve each do- 
main independently with communication done only at 
the boundaries. At the domain boundaries we compute 
the time derivatives of the characteristic fields in each 
domain. We then replace the time derivatives of the 
incoming characteristic fields at the boundary with the 
time derivatives of the outgoing characteristic fields of 
the neighboring domain. If there is no neighboring do- 
main at a particular boundary, the external boundary 
condition is applied as described in section 111 D . 



F. Solving elliptic equations 

In addition t o evo lving our hyperbolic system of evolu- 
tion equations (2.4), we may need to solve elliptic equa- 
tions in order to construct initial data for the Einstein- 
Klein-Gordon system or to enforce elliptic gauge condi- 
tions. Consider a linear elliptic equation of the form 



C{u{x)) = fix), 



(3.17) 



where u{x) is the solution we are seeking. This can be 
cast as a matrix problem where, unlike for FD, the matrix 
corresponding to the linear operator £ is full. In ID 
we solve this matrix equation directly, but for higher- 
dimensional problems, it will be more efficient to use an 
iterative method. 

A nonlinear elliptic equation such as the Hamiltonian 
constraint can be solved either by the methods described 
in ||38|| , or by linearizing the nonlinear system and iterat- 
ing the linearized equations until a solution is found. The 
latter method is employed in the work described here. 



G. Filtering 



E. Multiple domains 

In order to use a PSC method for problems of dimen- 
sion d greater than unity the computational domain must 
be sufiiciently simple that it can be mapped to I'' or 
X S*^ (where S*^ are two- spheres ) . For three dimen- 
sions, this typically means a cube, a sphere, or a spherical 
shell. If the computational domain is more complicated, 
then it must be decomposed into sub-domains that can 
each be mapped to one of these domains. For example, in 
two dimensions an L-shaped region can be decomposed 
into two adjacent rectangles. 

The binary black hole problem will need to be solved 
using multiple domains. Therefore we test our ability to 
handle multiple domains on our one-dimensional prob- 
lems. The use of multiple domains also provides a natu- 
ral way of making our code run in parallel. We use KeLP 



The errors in a spectral method are dominated by two 
types of terms of roughly equal magnitude. Truncation 
error arises from the neglect of the high-frequency terms 
that are not retained in the truncated series. Aliasing er- 
ror occurs because each neglected high-frequency mode 
is indistinguishable from some retained lower-frequency 
mode when sampled only at the collocation points; for ex- 
ample, the functions sin(7r2;/5) and sin(— 97ra;/5) take the 
same values on a grid of N points x G {0, 1, 2, . . . , TV— 1}. 
Because of aliasing error, power in the high-frequency 
mode, instead of being completely neglected, ends up 
contributing to the lower-frequency mode. 

When solving a nonlinear system of equations it be- 
comes important to control the aliasing error. This can 
be done by filtering the high-frequency modes of the re- 
tained series. For quadratic nonlinearities it is sufficient 
to zero the top third of the spectral coefficients to elim- 
inate aliasing Is^. In our ID evolutions, we have found 
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it necessary to filter only gauge variables that are com- 
puted from an elliptic equation. In effect we are smooth- 
ing the solutions to the elliptic gauge equations to elim- 
inate high-frequency noise. Our preliminary investiga- 
tions suggest that more extensive filtering may be re- 
quired to produce stable evolutions in 3D. 



IV. NUMERICAL RESULTS 

A. Schwarzschild black hole 

In this section we evolve a time-independent slicing 
of a Schwarzschild black hole. We begin our numeri- 
cal evolutions with initi al da ta corresponding to one of 
the slicings given in Sec. [ID. If the evolution equations 
are integrated exactly, the solution will remain time- 
independent. We can test the convergence of our method 
by measuring the deviation of the solution from the ini- 
tial data at a given coordinate time, or by measuring the 
constraint quantities (^.16D , which are zero for the ex- 
act solution. For all evolutions, the interior of the hole 
is excised, and no boundary condition is applied at the 
inner boundary because all characteristic fields are out- 
going (off the domain) there. In Table | we list the input 
parameters and the results for selected evolutions. 



1. Analytic gauge conditions 

The simplest gauge treatment is to fix the gauge vari- 
ables a and [3"^ to their initial values during the entire 
evolution. For example, one can begin the evolution with 
Kerr-Schild initial data ( 2.23) and set a a nd 13'^ according 
to the analytic expressions ( 3.23c-2.23d) for all time. In 
Figs, [l] and || we plot the norm of the Hamiltonian con- 
straint and the deviation of grr from the analytic solution 
versus time for such an evolution (run ^ from Table |) . 
Each plot shows results for several spatial resolutions Nr 
run at a fixed time resolution Ai — O.OOTAf . The features 
near t = lOM correspond to a small error pulse that be- 
gins at the outer boundary at t = 0, grows like as 
it propagates inwards, and eventually falls into the hole. 
After several crossing times, the evolution settles into a 
steady state that converges to the analytic solution as 
one increases the spatial resolution. We end the evolu- 
tions a,tt = llOOOM even though they clearly would have 
proceeded further. The convergence rate is exponential 
until machine roundoff errors dominate, as illustrated in 
Fig. ^. Repeating the evolutions shown in Figs. |^-^ for 
Painleve-GuUstrand initial data (run m from Table |) 
yields similar results. 

In Fig. 1^ we show the norm of the Hamiltonian con- 
straint for run ^ of Table |. This evolution is identical to 
run 1^ except the initial data, as well as the values of a and 
/3'" for all time, correspond to a time-independent har- 
monic slice of the Schwarzschild geometry ( 2.25| ). Rather 



TABLE L Input parameters for selected evolutions. For 
each evolution we list the initial data type, the locations r/M 
of the inner and outer boundaries, the outer boundary condi- 
tions on the incoming fields and on {(7^,17°}, the gauge 
conditions (including boundary conditions for elliptic equa- 
tions) , the time stepping algorithm, and the result of the evo- 
lution. 



Run ID'' 


IB 


OB 


BC"^ 


Lapse'^ 


Shift 


TS° Res* 




1 


KS 


1.9 


11.9 


F 


F 


C 


C 


R4 


Stb 




2 


KS 


1.9 


11.9 


c 


F 


C 


C 


R4 


Stb 




3 


KS 


1.9 


11.9 


F 


c 


C 


C 


R4 


Stb 




4 


KS 


1.9 


11.9 


c 


c 


C 


C 


R4 


Stb 




5 


KS 


1.75 


120 


C 


c 


C 


C 


R4 


Stb 




6 


KS 


1.9 


11.9 


F 


F 


S(cl)(cl) 


MD(X)(A) R4 


Exp 




7 


KS 


1.9 


11.9 


F 


F 


S(cl)(cl) 


MS(X)(A) 


R4 


Exp 




8 


KS 


1.9 


11.9 


c 


F 


S(cl)(cl) 


MS(X)(A) 


R4 


LG 




£ 


KS 


1.9 


11.9 


F 


c 


S(cl)(cl) 


MS(X)(A) 


R4 


Exp 




IC 


KS 


1.9 


11.9 


c 


c 


S(cl)(cl) 


MS(X)(A) 


R4 


LG 




n 


KS 


1.9 


11.9 


c 


c 


S(cl)(cl) 


MD(X)(A) R4 


LG 




12 


KS 


1.9 


11.9 


c 


c 


S(cl)(cl) 


MS(X)(A) 


BE 


LG 




13 


PG 


1.9 


11.9 


F 


F 


C 


c 


R4 


Stb 




14 


PG 


1.9 


11.9 


c 


F 


c 


c 


R4 


Stb 




15 


PG 


1.9 


11.9 


F 


c 


c 


c 


R4 


LG 




H 


PG 


1.9 


11.9 


c 


c 


c 


c 


R4 


LG 






PG 


1.75 


120 


C 


C 


c 


c 


R4 


Stb 




is 


PG 


1.9 


11.9 


F 


F 


S(c2)(F) 


MS(X)(R) 


R4 


Exp 




19 


PG 


1.9 


11.9 


c 


F 


S(c2)(F) 


MS(X)(R) 


R4 


Exp 




2C 


PG 


1.9 


11.9 


F 


c 


S(c2)(F) 


MS(X)(R) 


R4 


LG 




21 


PG 


1.9 


11.9 


c 


c 


S(c2)(F) 


MS(X)(R) 


R4 


LG 




22 


PG 


1.9 


11.9 


c 


c 


S(c2)(F) 


MS(X)(R) 


BE 


LG 




23 


H 


1.9 


11.9 


F 


F 


C 


C 


R4 


Exp 




"m 


H 


1.9 


3.9 


F 


F 


c 


c 


R4 


Stb 




25 


H 


1.9 


11.9 


C 


F 


c 


c 


R4 


Exp 




26 


H 


1.9 


11.9 


F 


C 


c 


c 


R4 


Exp 




27 


H 


1.9 


11.9 


C 


c 


c 


c 


R4 


QG 




28 


H 


1.9 


11.9 


C 


c 


c 


c 


BE 


QG 




29 


H 


1.75 


120 


C 


c 


c 


c 


R4 


LG 






H 


1.9 


11.9 


C 


c 


S(cl)(F) 


MS(X)(A) 


BE 


LG 






H 


1.9 


11.9 


C 


c 


S(ci)(F) 


MS(X)(A) 


R4 


LG 




32 


FH 


0.9 


10.9 


F 


F 


c 


c 


R4 


Exp 




33 


FH 


0.9 


10.9 


F 


C 


c 


c 


R4 


Exp 




34 


FH 


0.9 


10.9 


C 


F 


c 


c 


R4 


Exp 




35 


FH 


0.9 


10.9 


C 


C 


c 


c 


R4 


QG 




36 


FH 


0.9 


120 


C 


C 


c 


c 


R4 


LG 




37 


FH 


0.9 


6.9 


F 


F 


c 


c 


BE 


Stb 




38 


FH 


0.9 


7.9 


F 


F 


c 


c 


BE 


Exp 




39 


FH 


0.9 


10.9 


C 


C 


S(ci)(F) 


MS(X)(A) 


BE 


LG 



''PG: Painleve-GuUstrand; KS; Kerr-Schild; H; harmonic 
time; FH: fully harmonic 
'^F: freezing; C: constraint 

'^C: constant; S{ib){ob): stationary mean curvature; F: freeze; 
cn: c_ = — n at horizon or outer boundary 
"^C: constant; MS(i6)(o6): m inima l strain; MD{ib){ob): mini- 
mal distortion; X: equation (2.34); R: Robin; A: dtgr = 
'=R4 



4th-order Runge-Kutta; BE: backward Euler 
*Stb: stable; Exp: exponential growth; LG: linearly-growing 
gauge mode; QG: quadratically-growing gauge mode 
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FIG. 1. Long-term stability of the evolution of Kerr-Schild 
initial data, run |l| from Table l| Plotted is the £2 norm of the 
Hamiltonian constraint (2.16a) in units of M~'^ as a function 
of time for several spatial resolutions. The number of spectral 
coefficients Nr for each plot, starting at the top, is 12, 16, 20, 
24, 27, 32, 36, 40, 45, 48, 54, and 60. 
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FIG. 2. Norm of the error in grr as a function of time for 
the same evolutions shown in Fig. |l| 
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FIG. 3. Norms of the Hamiltonian constraint and errors 
in selected fundamental variables plotted as a function of the 
number of spectral coefficients A*',- a,t t — llOOOAf for the 
evolutions shown in Fig. hi. The quantities SKrr and Sfrrr are 
measured in units of A/- . The errors decrease exponentially 
with Nr. 



than settling to a steady state, the numerical solution 
grows exponentially at late times, eventually crashing 
the code. This is caused by a combination of high- 
frequency numerical instabilities, rapidly-growing gauge 
modes, and rapidly-growing constraint violating modes, 
all of which can be suppressed by appropriate changes in 
the evolution algorithm, as described below. Evolutions 
of fully harmonic initial data ( 2.26| ) behave similarly. It 
is not known why these instabilities are absent in evolu- 
tions of Kerr-Schild and Painleve-GuUstrand initial data. 
However, the dependence of stability on the choice of ini- 
tial data should not be too surprising if one thinks of the 
initial data as a background solution and the numerical 
evolution as a perturbation on this background: in gen- 
eral, modifying the background solution can change the 
stability of perturbations. 

For the evolutions shown in Fig. ^, freezing boundary 
conditions (2.2C) are imposed on the incoming charac- 



teristic fields U^, U~, and . One can suppress 
the constraint-violating modes seen in Fig. ^ by replac- 
ing the freezing boundary conditions on U^, U^, and U^^ 
with con straint boundary conditions as discussed in Sec- 
tion |T^. The resulting evolutions are shown in Figs. ^ 
and |6l Except for the evolution with Nr — 32 discussed 
below, the Hamiltonian constraint C settles to a steady 
state that converges exponentially to zero. The same is 
true for the other three constraints C^t, Crrr, and Cr- 
However, the metric quantities and other fundamental 
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FIG. 4. Norm of the Hamiltonian constraint (2.16a) as a 
function of time for several spatial resolutions for evolutions 
of harmonic initial data, run ^ of Table |. The number of 
spectral coefficients A''^ for each plot, starting at the top, is 
12, 15, 16, 18, 20, 24, 25, 27, 30, 32, and 36. 



variables grow approximately quadratically with time, 
eventually causing the simulations to terminate. Be- 
cause the constraints remain satisfied, we attribute this 
quadratic growth to a gauge mode. 

The Nr = 32 case shown in Figs. H and || suffers from 
high-frequency noise that grows exponentially in time. 
We have experimented with various methods of damp- 
ing this noise, including filtering the fundamental vari- 
ables after each time step and adding numerical dissipa- 
tion terms to the equations. However, we have obtained 
best results by changing our fourth order Runge-Kutta 
time-stepping algorithm to an implicit backwards Euler 
scheme, which is much more dissipative. Figures ^ and ^ 
show the results of this modification. The evolution now 
satisfies the constraints at late times for sufficiently fine 
resolution, but still suffers from a quadratically growing 
gauge mode that causes the coarser resolution runs to 
crash. This gauge mode can be suppressed by apply- 
ing active gauge conditions, as shown in Section IV A 2 
below. Evolutions of fully harmonic initial data j|2.2^ 
produce results similar to those shown in Figs. ^-|[ 

We note that even with analytic gauge conditions and 
freezing outer boundary conditions, evolutions of har- 
monic and fully harmonic initial data such as those shown 
in Fig. ^ become stable when the outer boundary is 
moved sufficiently close to the black hole (see runs 
and |3^ ). A similar dependence on the outer boundary 
location has also been reported by others 1^,^ • A pos- 
sible explanation for this is discussed briefly in |21|: For 




FIG. 5. Norm of the Hamiltonian constraint (2.16a) as a 
function of time for several spatial resolutions for evolutions 
of harmonic initial data, run ^ of Table ^. Constraint-based 
outer boundary conditions are imposed on (7°, Ut, and . 
Resolutions are the same as in Fig. 0. 
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FIG. 6. Error in grr versus time for the same evolutions 
shown in Fig. ^. The growth is quadratic in t at late times. 
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FIG. 7. Norm of the Hamiltonian constraint (2.16a) as a 
function of time for several spatial resolutions for evolutions 
of harmonic initial data, run of Table |. The evolutions 
are identical to those in Fig. |5| except a backwards Euler time 
stepping scheme is used. 
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FIG. 9. Norm of the Hamiltonian constraint (2.16a) as a 
function of time for several spatial resolutions Nr for evolu- 
tions of Painleve-GuUstrand initial data using elliptic gauge 
conditions, run ^ of Table |. The number of spectral coeffi- 
cients Nr for each plot, starting at the top, is 12, 16, 20, 24, 
27, 32, 36, 40, 45, and 48. 
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FIG. 8. Error in grr versus time for the same evolutions 
shown in Fig. |^. The growth is quadratic in t at late times. 



a nonzero shift vector, any unstable zero-speed modes 
present in the solution will propagate inward from the 
outer boundary with speed —/?'', growing as they propa- 
gate. If the domain is sufficiently small, these modes do 
not have time to grow appreciably before they are swal- 
lowed by the horizon. As discussed previously, we find 
that constraint boundary conditions suppress exponen- 
tially growing modes, and thus allow evolutions with a 
larger outer boundary radius (runs |^, |lj, 29 and |3^ ) . 



2. Elliptic gauge conditions 

Although choosing a time-independent f3^ and a is the 
simplest gauge condition to implement, for an evolving 
numerical solution such a choice does not actively enforce 
any particular coordinate condition. In fact, it is remark- 



able that many of the cases discussed in Section IV A 1 



remain stable when the coordinates experience small de- 
viations from the exact solution. For more than one black 
hole in three spatial dimensions, one will almost certainly 
need general gauge conditions designed to prevent large 
changes in the numerical solution of a stationary or quasi- 
stationary spacetime. 

Figure ^ shows the norm of the Hamiltonian constraint 
for an evolution of Painleve-GuUstrand initial data. The 
gauge variables (3^ and d are computed by solving the 
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FIG. 10. Error in g,.r versus time for the same evolutions 
shown in Fig. ^. For the highest resolution, the growth is only 
linear in t at late times. 



minimal strain and stationary mean curvature equa- 
tions ( |2.37| ) and (2.3!; ) after each time step. These elliptic 
equations require boundary conditions. We impose (2.34) 
and ( p.36| ) at the current location of the hori zon, which 
we recompute after every time step. For (2.36) we choose 
c_ = —2, which is the value of c_ at the horizon for the 
analytic solution ( ^.24 ). At the outer boundary, we set 
a = 1 and we impose a Robin condition ( 2.21 ) on with 
f^oo = 0, 71 = 1/2. As seen in the figure, the evolution 
remains stable and convergent. To achieve stability, we 
find it necessary to apply a simple 2/3 cutoff filter to a 
and each time they are computed, and to impose con- 
straint boundary conditions on and C/° (but not on 



C/^ ) . Figure ^ shows the error in grr for the same evo- 
lution. For the highest resolution, one can see a linearly- 
growing gauge mode. Although modes that grow linearly 
will eventually terminate a simulation, they pose no diffi- 
culty for long-term evolutions because a much longer run 
time can be achieved by a modest increase in resolution. 

Similar results for the case of harmonic initial data are 
shown in Figs, [ill and 12. The evolution is stable and con- 
vergent, and the rapidly- growing gauge mode that ter- 
minated the simulation in the case of time-independent 
f}^ and a (Section [VAl, Fig. ||) now grows only lin- 
early with time. As in the case shown in Fig. pi we 
use a backwards Euler scheme for time evolution. For a 
fourth-order Runge-Kutta time discretization, results are 
similar except the evolutions with Nr — 25, 30, and 32 
are unstable. 



10- 



IQi 



102 

t/M 



10= 



104 



FIG. 11. Norm of the Hamiltonian constraint (2.16a) as a 



function of time for several spatial resolutions for evolutions 
of harmonic initial data using elliptic gauge conditions, run ^ 
of Table |. Resolutions are the same as in Fig. 0. 
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FIG. 12. Error in grr versus time for the same evolutions 
shown in Fig. ^ The growth is only linear in t at late times. 
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B. Black hole plus scalar wave 



In this subsection, we add dynamics to our spherically 
symmetric spacetime by including a Klein- Gordon scalar 
field as a matter source, as described in Section [IE, To 



set up initial data, we first choose arbitrary initial values 
of n and $r- We then solve the Hamiltonian and mo- 
mentum constraints via the standard York-Lichnerowicz 
conformal decomposition |]5^ , using one of the time- 
independent representati ons o f the Schwarzschild geom- 
etry discussed in Section IIP as a background solution. 
Once we obtain the conformal factor and the trace-free 
longitudinal part of the extrinsic curvature from the con- 
straints, we reconstruct the fundamental variables. 

For the evolutions described here, the scalar field 11 
is initially a Gaussian centered at r = 20M with a 
width of 5M and an amplitude of 0.02/Af, and is 
zero everywhere. The outer boundary is located at 
r — 120M and the inner boundary is at 1.75M. Here 
M is the mass of the background solution, which is dif- 
ferent than the actual mass of the spacetime. We choose 
a Kerr-Schild background, analytic gauge conditions, and 
a fourth-order Runge-Kutta time stepping algorithm. At 
the outer boundary we apply freezing boundary condi- 
tions (2.20) to and , and constraint boundary 
conditions to [7°, [7° and Uj^ . There is no boundary 
condition imposed at the inner boundary. 

To demonstrate our ability to handle multiple domains, 
for this evolution we cover the entire domain with 8 equal- 
sized abutting subdomains, each using 45 spectral coef- 
ficients. At each domain boundary, the incoming char- 
acteristic quantities in each domain are set equal to the 
corresponding outgoing quantities of the neighboring do- 
main. 

In Fig. we plot the mass contained within radius r 
as a function of r for selected times. Initially the mass of 
the black hole is 0.97M and the mass of the entire space- 
time is 1.52M. The scalar field energy concentrated near 
r = 20M accounts for 0.55M. As the evolution proceeds, 
the initial Gaussian scalar field pulse divides into incom- 
ing and outgoing pieces. The outgoing piece propagates 
to infinity, while the incoming piece is partially refiected 
off the Schwarzschild potential and partially swallowed 
by the black hole. At t = 90M the mass of the black 
hole has reached its final value of 1.15M, and the ini- 
tial outgoing pulse and the reflected pulse have not yet 
reached the outer boundary of the domain. By t = 180M 
the remaining scalar radiation has left the domain. 

Figure ^ shows the coordinate radius and the areal ra- 
dius of the apparent horizon versus time. The area of the 
horizon increases between t = 15M and t — 30AI as the 
scalar field pulse falls into the black hole, and the areal 
radius asymptotically approaches the value 2.31M, which 
is twice the mass of the final black hole as expected. Un- 
like the areal radius, the coordinate radius decreases with 
time until t = lOM, after which it increases and even- 
tually asymptotes to r = 2.66Af . One can see from the 
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FIG. 13. Misner-Sharp mass, measured in units of AI, as 
a function of radius at selected times for an evolution of the 
Einstein-Klein-Gordon system. Here M is the mass of the 
background solution. See text for details. 
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FIG. 14. Coordinate radius (solid line) and areal radius 
(dotted line) of the apparent horizon as a function of time for 
the evolution shown in Fig. |l^. 
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FIG. 15. Norm of the Hamiltonian constraint (l2.16aD as a 
function of time for the evolution shown in Fig. |13| and for 
coarser evolutions with 12, 16, 18, 20, 24, 27, 30, 32, 36, and 
40 spectral coefficients per domain. 



figure that on the initial slice the coordinate r is nearly 
areal, as it would be for the Kerr-Schild background so- 
lution without a scalar field. At late times the deviation 
from an areal radial coordinate is large. 

Figure |l^ shows the norm of the Hamiltonian con- 
straint as a function of time for the same evolution, as 
well as for other evolutions with different spatial resolu- 
tions but the same At. The plots exhibit exponential con- 
vergence, even with nontrivial dynamics and multiple do- 
mains. There is, however, a small failure of convergence 
in the highest resolutions arou nd t — 120Af . This is be- 
cause the boundary condition ( 2.2C| ) , which is applied to 
U~ and at the outer boundary r;,, is not strictly cor- 
rect while a wave is passing through the boundary, and 
the error this introduces scales like r?,"^. We have ver- 
ified this by repeating the evolutions from Fig. |l5| with 
the inner boundary a factor of two closer, at 60M. This 
is shown in Fig, [iq. The resolution per subdomain is the 
same as in Fig. WK but we use 4 equal-sized subdomains 
instead of 8. The small nonconvergent feature in Fig. |l6| 
is approximately a factor of four larger than in Fig. 
and occurs at an earlier time, t — 60M, because the wave 
pulse reaches the outer boundary a factor of two earlier. 



V. DISCUSSION 

We have found that, at least for spherical symme- 
try, applying PSC methods to hyperbolic formulations of 



FIG. 16. Same as Fig. |l5| except each evolution has an 
outer boundary radius of 60Af instead of 120Af. 

general relativity can achieve stable evolutions of black 
holes with horizon excision. Excision itself is trivial as 
long as one uses a formulation in which all characteristic 
speeds are causal. Using realistic elliptic gauge condi- 
tions, our evolutions are limited only by linearly growing 
gauge modes that converge exponentially to zero with in- 
creasing resolution. These modes create no difficulty for 
long-term simulations because a small increase in resolu- 
tion enables one to run much farther in time. We note 
that even when errors grow exponentially in time (e.g.. 
Fig. U), the high accuracy provided by PSC allows us to 
evolve to times of hundreds or sometimes thousands of 
M. 

A hyperbolic formulation of Einstein's equations pro- 
vides a straightforward way in which to formulate and 
implement boundary conditions using the complete set 
of characteristic eigenfields provided by hyperbolicity. In 
principle, our method can be applied to any hyperbolic 
formulation, but so far we have only used the EC system. 
We have not investigated whether PSC can be used with 
non-hyperbolic formulations of Einstein's equations such 
as ADM, as this would require a different treatment of 
the boundary conditions. It is entirely possible, however, 
that one might find boundary conditions that result in 
stable evolutions for such a formulation. 

There has also been some concern about using hy- 
perbolic representations of general relativity with com- 
plicated gauge conditions. This is because hyperbolic 
formulations of Einstein's equations formally require the 
gauge quantities (shift and densitized lapse in the case of 
EC) to be prescribed functions of space and time, and not 
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FIG. 17. Two-dimensional illustration of multiple compu- 
tational domains that might be used to solve the binary black 
hole problem. Each hole is surrounded by a single domain in 
the shape of a spherical shell. Multiple cubical domains over- 
lap the spherical shells. There are sixteen shown here in two 
dimensions, with the cubes containing the black holes excised. 



evolved quantities that couple to the fundamental vari- 
ables. However, as long as the gauge variables are held 
fixed during each entire time step, as discussed by ||5^] 
and 1^^, we find no fundamental difficulty in applying 
elliptic gauge conditions during our simulations. 

By using the constraints as boundary conditions on 
the hyperbolic evolution equations, we have found that 
one can improve evolutions of the EC system. We have 
found similar improvement for finite-difference evolutions 
as well. Even in the general 3D case, applying constraint 
boundary conditions on the metric variables is straight- 
forward. However, casting the Hamiltonian and momen- 
tum constraints as boundary conditions may be more dif- 
ficult, because unlike in spherical symmetry, these con- 
straints involve contractions of derivatives of fundamen- 
tal variables. 

The numerical techniques discussed in this paper 
should be generalizable to three spatial dimensions. For 
PSC evolutions of two black holes with excised horizons it 
will be necessary to use multiple computational domains 
(see Fig. |l^). In spherical symmetry, we have shown how 
multiple domains can be easily implemented in a natural 
way by using characteristic fields to provide inter-domain 
boundary conditions. This method is directly applicable 
to abutting domains in 3D, and the extension from abut- 
ting domains to overlapping domains is straightforward 
pof . Work on 3D black hole evolutions using PSC is in 
progress. 



ACKNOWLEDGMENTS 

This work was supported in part by NSF grants PHY- 
9800737 and PHY-9900672 and NASA Grant NAG5-7264 
to Cornell University, and NSF grants PHY-9802571 and 
PHY-9988581 to Wake Forest University Computations 
were performed on the National Computational Science 
Alliance SGI Origin2000, and on the Wake Forest Univer- 



[9: 

[to; 

[11 
[12; 

[is: 

[14 

[is: 
[16: 

[ir 

[is: 

[19 
[20 
[21 

[22: 
[23: 

[24' 
[25^ 
[26: 

[27 

[28' 
[29: 



M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993). 

S. L. Shapiro, S. A. Teukolsky, and J. Winicour, Phys. 

Rev. D 52, 6982 (1995). 

S. F. P. Zwart and S. L. W. McMillan, Astrophys. J. 528, 
L17 (2000). 

E. E. Flanagan and S. A. Hughes, Phys. Rev. D 58, 4566 
(1998). 



W. Landry and S. A. Teukolsky, ^r-qc/9912004| (unpub- 
lished) . 

M. Shibata and K. Uryu, Phys. Rev. D 61, 064001 (2000). 
R. Arnowitt, S. Deser, and C. W. Misner, in Gravitation: 
An Introduction to Current Research, edited by L. Witten 
(Wiley, New York, 1962), pp. 227-265. 
S. FritteUi and O. Reula, Commun. Math. Phys. 166, 
221 (1994). 

Y. Choquet-Bruhat and J. W. York, Jr., C. R. Acad. Sci., 
Ser. I: Math. A321, 1089 (1995). 

A. Abrahams, A. Anderson, Y. Choquet-Bruhat, and 
J. W. York, Jr., Phys. Rev. Lett. 75, 3377 (1995). 
C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 
Lett. 75, 600 (1995). 

M. H. P. M. van Putten and D. M. Eardley, Phys. Rev. 
D 53, 3056 (1996). 

S. Frittelh and O. A. Reula, Phys. Rev. Lett. 76, 4667 
(1996). 

H. Friedrich, Class. Quantum Grav. 13, 1451 (1996). 

F. B. Estabrook, R. S. Robinson, and H. D. Wahlquist, 
Class. Quantum Grav. 14, 1237 (1997). 

A. Anderson, Y. Choquet-Bruhat, and J. W. York, Jr., 
Topol. Meth. Nonlin. Analy. 10, 353 (1997). 
M. Alcubierre, B. Brugmann, M. Miller, and W.-M. Suen, 
Phys. Rev. D 60, 064017 (1999). 

S. FritteUi and O. A. Reula, J. Math. Phys. 40, 5143 
(1999). 

A. Anderson and J. W. York, Jr., Phys. Rev. Lett. 82, 
4384 (1999). 

T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 
024007 (1999). 

M. A. Scheel et al, Phys. Rev. D 58, 044020 (1998). 
E. Seidel and W.-M. Suen, Phys. Rev. Lett. 69, 1845 
(1992). 

M. Alcubierre and B. F. Schutz, J. Comp. Phys. 112, 44 
(1994). 

M. A. Scheel et al., Phys. Rev. D 56, 6320 (1997). 

G. B. Cook et al, Phys. Rev. Lett. 80, 2512 (1998). 

C. Gundlach and P. Walker, Class. Quantum Grav. 16, 
991 (1999). 



L. Lehner, M. Huq, and D. Garrison, gr-qc/0004065 (un- 
published) . 

A. M. Abrahams et al, Phys. Rev. Lett. 80, 1812 (1998). 
N. T. Bishop, R. Gomez, L. Lehner, and J. Winicour, 
Phys. Rev. D 54, 6153 (1996). 



19 



R. Gomez et al., Phys. Rev. Lett. 80, 3915 (1998). 
N. T. Bishop et al, in Black Holes, Gravitational Ra- 
diation and the Universe, edited by B. R. Iyer and B. 
Bhawal (Kluwer, Dordrecht, 1998), Chap. 24. 
J. P. Boyd, Chebyshev and Fourier Spectral Methods 
( Springer- Ver lag, Berlin, Germany, 1989). 
C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. 
Zang, Spectral Methods in Fluid Dynamics (Springer- 
Verlag, Berlin, Germany, 1988). 

B. Fornberg, A Practical Guide to Pseudo spectral Meth- 
ods (Cambridge University Press, New York, 1996). 

S. Bonazzola, E. Gourgoulhon, and J. -A. Marck, J. Com- 
put. Appl. Math. 109, 892 (1999). 

J. Novak and J. M. Ibanez, Astrophys. J. 533, 392 (2000). 
S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. 
Rev. Lett. 82, 892 (1999). 

L. E. Kidder and L. S. Finn, ^r-qc/9911o"l^ (unpub- 
lished) . 

J. Frauendiener, J. Comput. Appl. Math. 109, 475 
(1999). 

C. Gundlach, Phys. Rev. D 57, 863 (1998). 

S. Bonazzola et al., in ICOSAHOM'95 Proceedings of the 
Third International Conference on Spectral and High Or- 
der Methods, edited by A. V. Ilin and L. R. Scott (Hous- 
ton Journal of Mathematics, Hous ton, 1996), pp. 3-19. 
R. A. Bartnik and A. H. Norton, |gr-qc/9904045| SIAM 
J. Sci. Comput. (to be published). 

R. L. Marsa and M. W. Choptuik, Phys. Rev. D 54, 4929 
(1996). 

P. Painleve, C. R. Acad. Sci. (Paris) 173, 677 (1921). 
A. GuUstrand, Arkiv. Mat Astron. Fys. 16, 1 (1922). 



gr-qc/0001069| (unpublished). 



K. Martel and E. Poisson, 

C. Bona and J. Masso, Phys. Rev. D 38, 2419 (1988). 
G. B. Cook and M. A. Scheel, Phys. Rev. D 56, 4775 
(1997). 

P. Anninos et al, Phys. Rev. D 51, 5562 (1995). 

L. L. Smarr and J. W. York, Jr., Phys. Rev. D 17, 2529 

(1978). 

D. Garfi nkle, C. Gund lach, J. Isenberg, and N. O. Mur- 



chadha, gr-qc/0003113 (unpublished). 

D. Garfinkle and C. Gundlach, Class. Quantum Grav. 

16, 4111 (1999). 

C. W. Misner, K. S. Thorne, and J. A. Wheeler, Grav- 
itation (W. H. Freeman and Company, New York, New 
York, 1973). 

S. B. Baden, D. Shalit, and R. B. Frost, KeLP User 



Guide, littp://www. cse.ucsd.edu/groups/hpcl/scg/kelp 



S. A. Orszag, J. Atmosph. Sci. 28, 1074 (1971). 
J. W. York, Jr., in Sources of Gravitational Radiation, 
edited by L. L. Smarr (Cambridge University Press, Cam- 
bridge, England, 1979), pp. 83-126. 
J. W. York, Jr., (private communication). 
C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. D 
56, 3405 (1997). 
[59] D. A. Kopriva, SIAM J. Sci. Stat. Comput. 10, 120 
(1989). 



20 



